clear
version 14.0
set more off
*net from http://www.stata-press.com/data/vgsg3/
*net install vgsg3
set scheme vg_s1m
capture log close
log using "..\logs\summary stats.txt", text replace

* Specify the year of cash rent data to use in the regression
local reg_year 2013
* Specify 1 if irrigated or 0 if nonirrigated
local practice 0

use "..\dataAnalysis\rent_climate_data", clear

if `practice'==0{
gen ln_cash_rent=ln(cash_rent_non`reg_year')
}
if `practice'==1{
gen ln_cash_rent=ln(cash_rent_irr`reg_year')
}

*------------------------------------------------------------------------
* Potential soils
*------------------------------------------------------------------------
local soils  slope om sar gypsum cec7 ec  ksat clay_perc silt_perc ///
			bulkDensity soc0_150 rootznemc ph_less6 ph_greater7dot5 hydgrp* text_*

local soils_nonlin  slope_* om_* sar_* gypsum_* cec7_* ec_*  ksat_* clay_perc_* silt_perc_* ///
			bulkDensity_* soc0_150_* rootznemc_* ph_less6 ph_greater7dot5 hydgrp* text_*

*------------------------------------------------------------------------
* Controls to include in the regression
*------------------------------------------------------------------------
local deficit water_deficit
local surplus water_surplus
local edd dday34C
local soils_controls soc0_150_4 bulkDensity_4 ec_2 ph_less6 ph_greater7dot5 slope_1     



*------------------------------------------------------------------------
* Summary Stats of Regression Variables
*------------------------------------------------------------------------
estpost tabstat cash_rent_non`reg_year' ln_cash_rent `deficit' `surplus' gdd `edd' `soils_controls' if cash_rent_non`reg_year'!=. & `deficit'!=., ///
		s(count mean sd min max) columns(statistics)
eststo summary_stats

esttab summary_stats using "..\tables\summary_stats.tex", ///
	 cells("mean(fmt(%9.2f)) sd(fmt(%9.2f))") noobs tex label nonumbers nomtitles replace ///
	 coeflabels(cash_rent_non`reg_year' "Cash Rent (\\$/acre)" ln_cash_rent "Log of Cash Rent" ///
	 `deficit' "Water Deficit (in)" `surplus' "Water Surplus (in)" ///
	 gdd "Hundreds of Growing Degree Days - GDD" dday34C "Extreme Degree Days - EDD"  ///
	 soc0_150_4 "Soil Organic Carbon (kg/m\$^2$) - SOC" bulkDensity_4 "Bulk Density" ///
	 ec_2 "Electrical conductivity - EC" ph_less6 "pH Less than 6" ph_greater7dot5 "pH Greater than 7.5" ///
	 slope_1 "Log of Slope")

corr water_deficit water_surplus `edd' ppt

corr gdd `edd' ET0

scatter gdd ET0
scatter `edd' ET0


*------------------------------------------------------------------------
* Summary Stats of Climate Projections
*------------------------------------------------------------------------
* Ensemble average change, sd across counties
estpost tabstat chg_water_deficit_rcp45AVG chg_water_surplus_rcp45AVG ///
				chg_gdd_rcp45AVG chg_dday34C_rcp45AVG ///
				, s(count mean sd ) columns(statistics)
eststo summary_projections

* sd across models
mean chg_water_deficit_rcp45AVG chg_water_surplus_rcp45AVG ///
				chg_gdd_rcp45AVG chg_dday34C_rcp45AVG 
eststo summary_projections

estout summary_projections, cells("mean sd b")

*------------------------------------------------------------------------
* Other Stats in the paper
*------------------------------------------------------------------------
corr ET0 vpd

summ chg_tmax_rcp45AVG chg_tmin_rcp45AVG
summ chg_tmax_rcp85AVG chg_tmin_rcp85AVG

summ chg_ppt_rcp45AVG chg_water_deficit_rcp45AVG, detail
*histogram chg_ppt_rcp45AVG
*histogram chg_water_deficit_rcp45AVG
			
*------------------------------------------------------------------------
* Graphs of climate projections across models
*------------------------------------------------------------------------
ren chg_water_deficit_* chg_deficit_*
ren chg_water_surplus_* chg_surplus_*
collapse (mean) chg_deficit_* chg_surplus_*	chg_gdd_* chg_dday34C_* 
gen obs=1

reshape long chg_deficit_rcp45 chg_surplus_rcp45 chg_gdd_rcp45 chg_dday34C_rcp45 ///
		chg_deficit_rcp85 chg_surplus_rcp85 chg_gdd_rcp85 chg_dday34C_rcp85, i(obs) j(model) string

gen marker_label=""
replace marker_label="A" if model=="AVG"	

label variable chg_deficit_rcp45 "Change Water Deficit" 
label variable chg_surplus_rcp45 "Change Water Surplus" 
label variable chg_gdd_rcp45 "Change GDD" 
label variable chg_dday34C_rcp45 "Change EDD" 	
label variable chg_deficit_rcp85 "Change Water Deficit" 
label variable chg_surplus_rcp85 "Change Water Surplus" 
label variable chg_gdd_rcp85 "Change GDD" 
label variable chg_dday34C_rcp85 "Change EDD" 	
		
graph matrix chg_deficit_rcp45 chg_surplus_rcp45 chg_gdd_rcp45 chg_dday34C_rcp45, ///
	 half msymbol(Oh) mcolor(gs4) msize(medlarge) ///
	 mlabel(marker_label) mlabp(0) mlabc(cranberry) iscale(1.25)

graph export "..\figures\distribution_projections_rcp45.pdf", replace	


summ *rcp45 if model=="AVG" 
summ *rcp45 

graph matrix chg_deficit_rcp85 chg_surplus_rcp85 chg_gdd_rcp85 chg_dday34C_rcp85, ///
	 half msymbol(Oh) mcolor(gs4) msize(medlarge) ///
	 mlabel(marker_label) mlabp(0) mlabc(cranberry)  iscale(1.25)

graph export "..\figures\distribution_projections_rcp85.pdf", replace

*------------------------------------------------------------------------
* Summary Stats of Climate Projections
* I don't end up including these tables in the paper, but they are
* still useful to look at
*------------------------------------------------------------------------
use "..\dataAnalysis\rent_climate_data", clear
if `practice'==0{
gen ln_cash_rent=ln(cash_rent_non`reg_year')
}
if `practice'==1{
gen ln_cash_rent=ln(cash_rent_irr`reg_year')
}

ren chg_water_deficit_* chg_deficit_*
ren chg_water_surplus_* chg_surplus_*
collapse (mean) chg_deficit_* chg_surplus_*	chg_gdd_* chg_dday34C_* ///
		(sd) sd_deficit_rcp45AVG=chg_deficit_rcp45AVG sd_surplus_rcp45AVG=chg_surplus_rcp45AVG ///
		sd_gdd_rcp45AVG=chg_gdd_rcp45AVG sd_dday34C_rcp45AVG=chg_dday34C_rcp45AVG ///
		sd_deficit_rcp85AVG=chg_deficit_rcp85AVG sd_surplus_rcp85AVG=chg_surplus_rcp85AVG ///
		sd_gdd_rcp85AVG=chg_gdd_rcp85AVG sd_dday34C_rcp85AVG=chg_dday34C_rcp85AVG
		
local future_scenarios "rcp45 rcp85"
		
foreach fs of local future_scenarios {		
* Average change		
mean chg_deficit_`fs'AVG chg_surplus_`fs'AVG ///
				chg_gdd_`fs'AVG chg_dday34C_`fs'AVG
eststo mean_chg_`fs'
drop chg_deficit_`fs'AVG chg_surplus_`fs'AVG ///
				chg_gdd_`fs'AVG chg_dday34C_`fs'AVG

* Standard Deviation across counties
ren sd_deficit_`fs'AVG chg_deficit_`fs'AVG
ren sd_surplus_`fs'AVG chg_surplus_`fs'AVG
ren sd_gdd_`fs'AVG chg_gdd_`fs'AVG
ren sd_dday34C_`fs'AVG	chg_dday34C_`fs'AVG		
mean chg_deficit_`fs'AVG chg_surplus_`fs'AVG ///
				chg_gdd_`fs'AVG chg_dday34C_`fs'AVG			
eststo sd_counties_`fs'

* Standard Deviation across models
drop chg_deficit_`fs'AVG chg_surplus_`fs'AVG ///
				chg_gdd_`fs'AVG chg_dday34C_`fs'AVG
egen chg_deficit_`fs'AVG=rowsd(chg_deficit_`fs'*)	
egen chg_surplus_`fs'AVG=rowsd(chg_surplus_`fs'*)	
egen chg_gdd_`fs'AVG=rowsd(chg_gdd_`fs'*)	
egen chg_dday34C_`fs'AVG=rowsd(chg_dday34C_`fs'*)	
mean chg_deficit_`fs'AVG chg_surplus_`fs'AVG ///
				chg_gdd_`fs'AVG chg_dday34C_`fs'AVG			
eststo sd_models_`fs'	

* Minimum across models
drop chg_deficit_`fs'AVG chg_surplus_`fs'AVG ///
				chg_gdd_`fs'AVG chg_dday34C_`fs'AVG
egen chg_deficit_`fs'AVG=rowmin(chg_deficit_`fs'*)	
egen chg_surplus_`fs'AVG=rowmin(chg_surplus_`fs'*)	
egen chg_gdd_`fs'AVG=rowmin(chg_gdd_`fs'*)	
egen chg_dday34C_`fs'AVG=rowmin(chg_dday34C_`fs'*)	
mean chg_deficit_`fs'AVG chg_surplus_`fs'AVG ///
				chg_gdd_`fs'AVG chg_dday34C_`fs'AVG			
eststo min_models_`fs'

* Maximum across models
drop chg_deficit_`fs'AVG chg_surplus_`fs'AVG ///
				chg_gdd_`fs'AVG chg_dday34C_`fs'AVG
egen chg_deficit_`fs'AVG=rowmax(chg_deficit_`fs'*)	
egen chg_surplus_`fs'AVG=rowmax(chg_surplus_`fs'*)	
egen chg_gdd_`fs'AVG=rowmax(chg_gdd_`fs'*)	
egen chg_dday34C_`fs'AVG=rowmax(chg_dday34C_`fs'*)	
mean chg_deficit_`fs'AVG chg_surplus_`fs'AVG ///
				chg_gdd_`fs'AVG chg_dday34C_`fs'AVG			
eststo max_models_`fs'

esttab mean_chg_`fs' sd_models_`fs' min_models_`fs' max_models_`fs' ///
		using "..\tables\summary_projections_`fs'.tex", b(%9.2f) ///
		 not nostar noobs tex label nonumbers replace ///
		 mtitles("Ensemble Average" "SD" "Min" "Max") ///
		 coeflabels(chg_deficit_`fs'AVG "Change in Water Deficit" ///
		 chg_surplus_`fs'AVG "Change in Water Surplus" ///
		 chg_gdd_`fs'AVG "Change in GDD" ///
		 chg_dday34C_`fs'AVG "Change in EDD")	

}

		 
log close
